function pdd (L)
 
   ## usage: pdd (L)
   ## 
   ## process predictands (stations)

   global PAR PRJ ids vars names lons lats alts undefs clim

   ldata = fullfile("data", PRJ, PAR.pdd) ;
   METAFILE = fullfile(ldata, "meta.txt") ;
   ofile = fullfile(ldata, "pdd.ob") ;
   if (exist(ofile, "file") && !L), return ; endif

   f = dbstack ; printf("\n\t###   %s   ###\n\n", f(1).name) ;

   pkg load general miscellaneous nan statistics

   [ids, vars, names, lons, lats, alts, undefs, J] = meta(METAFILE) ;
   VAR = sunique(vars) ;

   pdd = [] ;
   for j = 1:length(VAR)

      if !any(J{j}), continue ; endif 

      v = VAR{j} ;

      s = do_read(j, v, J{j}) ;
      s = fill_date(s) ;	#  fix broken observations
      s = do_norm(s) ;

      pdd = join(pdd, s) ;

   endfor

   [pdd.nr pdd.nc] = size(pdd.x) ;
   clim.vars = pdd.vars ; clim.ids = pdd.ids ; clim.names = pdd.names ; clim.lons = pdd.lons ; clim.lats = pdd.lats ;
   clim.isP = pdd.isP ; clim.ispos = pdd.ispos ;

   save(cfile = fullfile(ldata, "clim.ob"), "clim") ;
   printf("pdd: --> %s\n", cfile) ;

   save(ofile, "pdd") ;
   printf("pdd: --> %s\n", ofile) ;
   
endfunction


function s = do_read (j, v, J)

   ## usage: s = do_read (j, v, J)
   ## 
   ## 

   global ids vars names lons lats alts undefs
   global PAR PRJ

   ldata = fullfile("data", PRJ, PAR.pdd) ;

   Jv = cellfun(@(x) ~isempty(x), strfind(vars, v)) ; vn = tolower(names(Jv)) ;

   file = fullfile(ldata, [v ".csv"]) ;
   printf("pdd: <-- %s\n", file) ;
   s = read_stat(file) ;

   isP = any(cellfun(@(x) !isempty(x), strfind(vn, "precip"))) ;
   ispos = isP | any(cellfun(@(x) !isempty(x), strfind(vn, "radiation"))) ;
   s.isP = repmat(isP, 1, columns(s.x)) ;
   s.ispos = repmat(ispos, 1, columns(s.x)) ;

   v = strrep(v, ".", "_") ;
   s.vars = repmat({v}, 1, columns(s.x)) ;

   [JP JJ] = pick_meta(ids, vars, s.ids, v) ;

   if !all(strcmp(s.ids'(JJ), ids(JP)(JJ)))
      warning off backtrace
      warning("xds:xds", "pdd: unequal id ordering between %s and %s", file, fullfile(ldata, "meta.txt")) ;
   endif

   s.ids = ids(JP)(JJ)' ;
   s.vars = vars(JP)(JJ)' ;
   s.names = names(JP)(JJ)' ;
   s.lons = str2double(lons(JP)(JJ)') ;
   s.lats = str2double(lats(JP)(JJ)') ;
   s.alts = str2double(alts(JP)(JJ)') ;
   s.udf = str2double(undefs(JP)(JJ)') ;

   s.x = s.x(:,JJ) ; s.x = lim2nan(s.x, s.udf) ;
   s.isP = s.isP(:,JJ) ;
   s.ispos = s.ispos(:,JJ) ;

   datagaps(s) ;

endfunction


function datagaps (s)

   ## usage:  res = datagaps (s)
   ##
   ## check overlap of data

   global BATCH PAR PRJ

   ldata = fullfile("data", PRJ, PAR.pdd) ;

   if strcmp(PAR.optptr, "manual") || true # FIXME: needed by skl in cal
      sw = selper(s, sdate(s.id, PAR.xdscal(1,:), PAR.xdscal(2,:))) ;
   else
      sw = selper(s, sdate(s.id, PAR.cper(1,:), PAR.cper(2,:))) ;
   endif
   Idfd = !isnan(sw.x) ;
   I = Idfd'* Idfd ;

   if any(II = (I(:) < PAR.Ccases))
      disp(s.ids) ; disp(I) ;
      h = disp_dfd(sw) ;
      if BATCH
	 pause(10) ;
	 printf("dsc: --> %s\n", pfile = fullfile(ldata, [s.vars{1} ".dfd.png"])) ;
	 print(pfile) ;
	 close(h) ;
      endif
      error("pdd : not enough overlapping calibration data (%d < %d).", I(:)(II)(1), PAR.Ccases)
   endif 

   if BATCH && PAR.dbg
      h = disp_dfd(sw) ;
      printf("dsc: --> %s\n", pfile = fullfile(ldata, [s.vars{1} ".dfd.png"])) ;
      print(pfile) ;
      close(h) ;      
   endif

endfunction


function v = do_norm (u)

   ## usage: v = do_norm (u)
   ## 
   ## 

   global PAR clim

   v = u ;

   vn = v.vars{1} ;

   v = pre_pare(v) ;

   ##vw = selper(v, sdate(v.id, PAR.cper(1,:), PAR.cper(2,:))) ;
   vw = v ;

   if !PAR.anc
      v.a = v.x ;
      vw.a = vw.x ;
      [vw.z, clim.(vn).sm, clim.(vn).ss] = anom(vw.a) ;
      v.z = anom(v.a, clim.(vn).sm, clim.(vn).ss) ;
   else
      if all(v.ispos)
	 if PAR.dbg, printf(["pdd: " vn " treated as positive.\n"]) ; endif
	 [clim.(vn).ps, clim.(vn).pp, vw] = anncyclog(vw) ;
	 v = anncyclog(v, clim.(vn).ps, clim.(vn).pp) ;
	 v.z = pfy(v.z, clim.(vn).pp(1,:)) ;
      else
	 [clim.(vn).ps vw] = anncyc(vw) ;
	 [clim.(vn).pp clim.(vn).ppinv] = prbfit(vw.a) ;
	 v = anncyc(v, clim.(vn).ps) ;
	 v.z = prbfit(v.a, clim.(vn).pp) ;
      endif 
   endif

   v.x = u.x ;			# restore x

endfunction
